package computation;

import java.util.ArrayList;
import java.util.HashSet;

import triangulation.DelaunayMapping;
import util.GeoEdge;
import util.GeoTriangle;
import util.GeoVertex;
import util.Matrix_3X3;

public class Curvatures
{
	public static double computeMeanCurvature(ArrayList<GeoEdge> edges, DelaunayMapping mapping, int smoothness)
	{
		double curvature = 0.0;
		for (GeoEdge e: edges)
		{
			GeoTriangle[] incident = new GeoTriangle[2];
			incident[0] = null;
			incident[1] = null;
			
			HashSet<Integer> startTriangles = e.start.trianglekeys;
			HashSet<Integer> endTriangles = e.end.trianglekeys;
			
			for (Integer stk: startTriangles)
			{
				if (endTriangles.contains(stk))
				{
					if (incident[0] == null) incident[0] = mapping.triangles[stk];
					else if (incident[1] == null) incident[1] = mapping.triangles[stk];
					else break;
				}
			}
			
			double total_angle = 0.0; // B(e) from Steiner's paper
			// TODO: only deal with edges that have two triangles for now, find out how to handle single triangle cases later
			if (incident[0] != null && incident[1] != null)
			{
				// compute the dot product
				GeoVertex g1_0 = mapping.getVertex(incident[0].vertexkeys[0], smoothness);
				GeoVertex g1_1 = mapping.getVertex(incident[0].vertexkeys[1], smoothness);
				GeoVertex g1_2 = mapping.getVertex(incident[0].vertexkeys[2], smoothness);
				GeoVertex g2_0 = mapping.getVertex(incident[1].vertexkeys[0], smoothness);
				GeoVertex g2_1 = mapping.getVertex(incident[1].vertexkeys[1], smoothness);
				GeoVertex g2_2 = mapping.getVertex(incident[1].vertexkeys[2], smoothness);
				float[] n1 = Geometry.computeSurfaceNormal(g1_0, g1_1, g1_2);
				float[] n2 = Geometry.computeSurfaceNormal(g2_0, g2_1, g2_2);
				float dotproduct = 0.0f;
				for(int i = 0; i < 3; i++) dotproduct += (n1[i] * n2[i]);
				double angle = Math.acos(dotproduct / (Geometry.vectorMagnitude(n1) * Geometry.vectorMagnitude(n2)));
//				total_angle += (angle/2);
				total_angle += Math.toDegrees(angle)/2;
			}
			curvature += (e.length()*total_angle);
		}
		return curvature/edges.size();
	}
	
	public static double[] computeAnisotropicCurvature(ArrayList<GeoEdge> edges, DelaunayMapping mapping, int smoothness)
	{
		double[] eigenvalues = new double[3];
		double[][] curvature_tensor = new double[3][3];
		for (GeoEdge e: edges)
		{
			GeoTriangle[] incident = new GeoTriangle[2];
			incident[0] = null;
			incident[1] = null;
			
			HashSet<Integer> startTriangles = e.start.trianglekeys;
			HashSet<Integer> endTriangles = e.end.trianglekeys;
			for (Integer stk: startTriangles)
			{
				if (endTriangles.contains(stk))
				{
					if (incident[0] == null) incident[0] = mapping.triangles[stk];
					else if (incident[1] == null) incident[1] = mapping.triangles[stk];
					else break;
				}
			}

			// TODO: only deal with edges that have two triangles for now, find out how to handle single triangle cases later
			if (incident[0] != null && incident[1] != null)
			{				
				// now we need to compute dir (tensor) dir^t
				double[] dir = e.getDirectionVector_Normalized();
				
				// compute B(e) from steiners paper
				GeoVertex g1_0 = mapping.getVertex(incident[0].vertexkeys[0], smoothness);
				GeoVertex g1_1 = mapping.getVertex(incident[0].vertexkeys[1], smoothness);
				GeoVertex g1_2 = mapping.getVertex(incident[0].vertexkeys[2], smoothness);
				
				GeoVertex g2_0 = mapping.getVertex(incident[1].vertexkeys[0], smoothness);
				GeoVertex g2_1 = mapping.getVertex(incident[1].vertexkeys[1], smoothness);
				GeoVertex g2_2 = mapping.getVertex(incident[1].vertexkeys[2], smoothness);
				
				float[] n1 = Geometry.computeSurfaceNormal(g1_0, g1_1, g1_2);
				float[] n2 = Geometry.computeSurfaceNormal(g2_0, g2_1, g2_2);
				float dotproduct = 0.0f;
				for(int i = 0; i < 3; i++) dotproduct += (n1[i] * n2[i]);
				double angle = Math.acos(dotproduct / (Geometry.vectorMagnitude(n1) * Geometry.vectorMagnitude(n2)));
				double total_angle = Math.toDegrees(angle);
				
				// now we compute the sign of B(e) by taking the normals of t1 and t2, and the two points (p1, p2) not connected by the edge e
				// convex edge = positive sign, concave edge = negative sign
				int sign = 1;
				GeoVertex start = e.start;
				GeoVertex end = e.end;
				if (e.isParametric()) end = e.extension; // need the vertex of the end of the full edge in order to determine p1 & p2

				// p1 & p2 are the two vertices of t1 & t2 which do not line on the shared edge
				GeoVertex p1 = null;
				GeoVertex p2 = null;
				if ((start.id == g1_0.id && end.id == g1_1.id) || (start.id == g1_1.id && end.id == g1_0.id)) p1 = g1_2;
				if ((start.id == g1_0.id && end.id == g1_2.id) || (start.id == g1_2.id && end.id == g1_0.id)) p1 = g1_1;
				if ((start.id == g1_2.id && end.id == g1_1.id) || (start.id == g1_1.id && end.id == g1_2.id)) p1 = g1_0;
				
				if ((start.id == g2_0.id && end.id == g2_1.id) || (start.id == g2_1.id && end.id == g2_0.id)) p2 = g2_2;
				if ((start.id == g2_0.id && end.id == g2_2.id) || (start.id == g2_2.id && end.id == g2_0.id)) p2 = g2_1;
				if ((start.id == g2_2.id && end.id == g2_1.id) || (start.id == g2_1.id && end.id == g2_2.id)) p2 = g2_0;
				
				// find out if our edge is concave or convex
				double[] diff = Geometry.vectorMinus(p2, p1);
				float dot = 0.0f;
				for(int i = 0; i < 3; i++) dot += (diff[i] * n1[i]);
				
				if (dot > 0.0) sign = -1; // concave
				
				double factor = e.length()*total_angle*sign;
				
				curvature_tensor[0][0] += (dir[0]*dir[0]*factor);
				curvature_tensor[0][1] += (dir[0]*dir[1]*factor);
				curvature_tensor[0][2] += (dir[0]*dir[2]*factor);
				
				curvature_tensor[1][0] += (dir[1]*dir[0]*factor);
				curvature_tensor[1][1] += (dir[1]*dir[1]*factor);
				curvature_tensor[1][2] += (dir[1]*dir[2]*factor);
				
				curvature_tensor[2][0] += (dir[2]*dir[0]*factor);
				curvature_tensor[2][1] += (dir[2]*dir[1]*factor);
				curvature_tensor[2][2] += (dir[2]*dir[2]*factor);
			}
		}
		
		eigenvalues = Matrix_3X3.computeEigenvalues(curvature_tensor);
		return eigenvalues;
	}
	
}
